home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / ASTRNOMY / AA_51.ZIP / PRECESS.C < prev    next >
C/C++ Source or Header  |  1993-02-09  |  6KB  |  239 lines
  1. /* Precession of the equinox and ecliptic
  2.  * from epoch Julian date J to or from J2000.0
  3.  *
  4.  * Program by Steve Moshier.
  5.  *
  6.  *
  7.  * IAU Coefficients are from:
  8.  * J. H. Lieske, T. Lederle, W. Fricke, and B. Morando,
  9.  * "Expressions for the Precession Quantities Based upon the IAU
  10.  * (1976) System of Astronomical Constants,"  Astronomy and
  11.  * Astrophysics 58, 1-16 (1977).
  12.  *
  13.  * Newer formulas that cover a much longer time span are from:
  14.  * J. Laskar, "Secular terms of classical planetary theories
  15.  * using the results of general theory," Astronomy and Astrophysics
  16.  * 157, 59070 (1986).
  17.  *
  18.  * See also:
  19.  * P. Bretagnon and G. Francou, "Planetary theories in rectangular
  20.  * and spherical variables. VSOP87 solutions," Astronomy and
  21.  * Astrophysics 202, 309-315 (1988).
  22.  *
  23.  * Laskar's expansions are said by Bretagnon and Francou
  24.  * to have "a precision of about 1" over 10000 years before
  25.  * and after J2000.0 in so far as the precession constants p^0_A
  26.  * and epsilon^0_A are perfectly known."
  27.  *
  28.  * Bretagnon and Francou's expansions for the node and inclination
  29.  * of the ecliptic were derived from Laskar's data but were truncated
  30.  * after the term in T**6. I have recomputed these expansions from
  31.  * Laskar's data, retaining powers up to T**10 in the result.
  32.  *
  33.  * The following table indicates the differences between the result
  34.  * of the IAU formula and Laskar's formula using four different test
  35.  * vectors, checking at J2000 plus and minus the indicated number
  36.  * of years.
  37.  *
  38.  *   Years       Arc
  39.  * from J2000  Seconds
  40.  * ----------  -------
  41.  *        0      0
  42.  *      100    .006    
  43.  *      200     .006
  44.  *      500     .015
  45.  *     1000     .28
  46.  *     2000    6.4
  47.  *     3000   38.
  48.  *    10000 9400.
  49.  */
  50.  
  51.  
  52. /* Precession coefficients taken from Laskar's paper: */
  53. static double pAcof[] = {
  54.  -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
  55.  -0.235316, 0.07732, 111.1971, 50290.966 };
  56.  
  57. /* Node and inclination of the earth's orbit computed from
  58.  * Laskar's data as done in Bretagnon and Francou's paper:
  59.  */
  60. static double nodecof[] = {
  61. 6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 6.3190131e-10, 
  62. -3.48388152e-9, -1.813065896e-7, 2.75036225e-8, 7.4394531426e-5,
  63. -0.042078604317, 3.052112654975 };
  64.  
  65. static double inclcof[] = {
  66. 1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11, 
  67. -5.4000441e-11, 1.32115526e-9, -5.998737027e-7, -1.6242797091e-5,
  68.  0.002278495537, 0.0 };
  69.  
  70. extern double J2000; /* = 2451545.0, 2000 January 1.5 */
  71. extern double STR; /* = 4.8481368110953599359e-6 radians per arc second */
  72. extern double coseps, sineps; /* see epsiln.c */
  73.  
  74. /* Subroutine arguments:
  75.  *
  76.  * R = rectangular equatorial coordinate vector to be precessed.
  77.  *     The result is written back into the input vector.
  78.  * J = Julian date
  79.  * direction =
  80.  *      Precess from J to J2000: direction = 1
  81.  *      Precess from J2000 to J: direction = -1
  82.  * Note that if you want to precess from J1 to J2, you would
  83.  * first go from J1 to J2000, then call the program again
  84.  * to go from J2000 to J2.
  85.  */
  86. extern int epsiln();
  87.  
  88. int precess( R, J, direction )
  89. double R[], J;
  90. int direction;
  91. {
  92. double sinth, costh, sinZ, cosZ, sinz, cosz;
  93. double A, B, T, Z, z, TH, pA, W;
  94. double x[3];
  95. double *p;
  96. int i;
  97. double sin(), cos(), fabs();
  98.  
  99. if( J == J2000 )
  100.     return(0);
  101. /* Each precession angle is specified by a polynomial in
  102.  * T = Julian centuries from J2000.0.  See AA page B18.
  103.  */
  104. T = (J - J2000)/36525.0;
  105.  
  106. if( fabs(T) > 2.0 )
  107.     goto laskar;
  108.  
  109. Z =  (( 0.017998*T + 0.30188)*T + 2306.2181)*T*STR;
  110. z =  (( 0.018203*T + 1.09468)*T + 2306.2181)*T*STR;
  111. TH = ((-0.041833*T - 0.42665)*T + 2004.3109)*T*STR;
  112.  
  113. sinth = sin(TH);
  114. costh = cos(TH);
  115. sinZ = sin(Z);
  116. cosZ = cos(Z);
  117. sinz = sin(z);
  118. cosz = cos(z);
  119. A = cosZ*costh;
  120. B = sinZ*costh;
  121.  
  122. if( direction < 0 )
  123.     { /* From J2000.0 to J */
  124.     x[0] =    (A*cosz - sinZ*sinz)*R[0]
  125.             - (B*cosz + cosZ*sinz)*R[1]
  126.                       - sinth*cosz*R[2];
  127.  
  128.     x[1] =    (A*sinz + sinZ*cosz)*R[0]
  129.             - (B*sinz - cosZ*cosz)*R[1]
  130.                       - sinth*sinz*R[2];
  131.  
  132.     x[2] =              cosZ*sinth*R[0]
  133.                       - sinZ*sinth*R[1]
  134.                            + costh*R[2];
  135.     }
  136. else
  137.     { /* From J to J2000.0 */
  138.     x[0] =    (A*cosz - sinZ*sinz)*R[0]
  139.             + (A*sinz + sinZ*cosz)*R[1]
  140.                       + cosZ*sinth*R[2];
  141.  
  142.     x[1] =   -(B*cosz + cosZ*sinz)*R[0]
  143.             - (B*sinz - cosZ*cosz)*R[1]
  144.                       - sinZ*sinth*R[2];
  145.  
  146.     x[2] =             -sinth*cosz*R[0]
  147.                       - sinth*sinz*R[1]
  148.                            + costh*R[2];
  149.     }    
  150. goto done;
  151.  
  152. laskar:
  153.  
  154. /* Implementation by elementary rotations using Laskar's expansions.
  155.  * First rotate about the x axis from the initial equator
  156.  * to the ecliptic. (The input is equatorial.)
  157.  */
  158. if( direction == 1 )
  159.     epsiln( J ); /* To J2000 */
  160. else
  161.     epsiln( J2000 ); /* From J2000 */
  162. x[0] = R[0];
  163. z = coseps*R[1] + sineps*R[2];
  164. x[2] = -sineps*R[1] + coseps*R[2];
  165. x[1] = z;
  166.  
  167. /* Precession in longitude
  168.  */
  169. T /= 10.0; /* thousands of years */
  170. p = pAcof;
  171. pA = *p++;
  172. for( i=0; i<9; i++ )
  173.     pA = pA * T + *p++;
  174. pA *= STR * T;
  175.  
  176. /* Node of the moving ecliptic on the J2000 ecliptic.
  177.  */
  178. p = nodecof;
  179. W = *p++;
  180. for( i=0; i<10; i++ )
  181.     W = W * T + *p++;
  182.  
  183. /* Rotate about z axis to the node.
  184.  */
  185. if( direction == 1 )
  186.     z = W + pA;
  187. else
  188.     z = W;
  189. B = cos(z);
  190. A = sin(z);
  191. z = B * x[0] + A * x[1];
  192. x[1] = -A * x[0] + B * x[1];
  193. x[0] = z;
  194.  
  195. /* Rotate about new x axis by the inclination of the moving
  196.  * ecliptic on the J2000 ecliptic.
  197.  */
  198. p = inclcof;
  199. z = *p++;
  200. for( i=0; i<10; i++ )
  201.     z = z * T + *p++;
  202. if( direction == 1 )
  203.     z = -z;
  204. B = cos(z);
  205. A = sin(z);
  206. z = B * x[1] + A * x[2];
  207. x[2] = -A * x[1] + B * x[2];
  208. x[1] = z;
  209.  
  210. /* Rotate about new z axis back from the node.
  211.  */
  212. if( direction == 1 )
  213.     z = -W;
  214. else
  215.     z = -W - pA;
  216. B = cos(z);
  217. A = sin(z);
  218. z = B * x[0] + A * x[1];
  219. x[1] = -A * x[0] + B * x[1];
  220. x[0] = z;
  221.  
  222. /* Rotate about x axis to final equator.
  223.  */
  224. if( direction == 1 )
  225.     epsiln( J2000 );
  226. else
  227.     epsiln( J );
  228. z = coseps * x[1] - sineps * x[2];
  229. x[2] = sineps * x[1] + coseps * x[2];
  230. x[1] = z;
  231.  
  232.  
  233. done:
  234.  
  235. for( i=0; i<3; i++ )
  236.     R[i] = x[i];
  237. return(0);
  238. }
  239.